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We investigate the electronic structure induced by wedge-disclinations (conical singularities) in a 
honeycomb lattice model realizing Chern numbers 7 = ±1. We establish a correspondence between 
the bound state of (i) an isolated $o/2-flux, (ii) an isolated pentagon (n = 1) or heptagon (n = — 1) 
defect with an external flux of magnitude n7$o/4 through the center and (iii) an isolated square 
or octagon defect without external flux, where $0 = h/e is the flux quantum. Due to the above 
correspondence, the existence of isolated electronic states bound to the disclinations is robust against 
various perturbations. These results are also generalized to graphene-based time-reversal invariant 
topological insulators. 
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When a magnetic flux tube is inserted into a three- 
dimensional time-reversal invariant topological insulator 
(TI) [1—3], a single Kramer's pair of counter propagat- 
ing gapless modes along the tube appears if the flux is a 
half-integer multiple of the flux quantum $0 = h/e — 2tt 
[4]. Physical manifestations of this "wormhole effect" 
[5] include conductance oscillations with period $0 in TI 
nanowires [6-8] as well as protected gapless modes along 
certain dislocations [9]. In the latter case, the perfectly 
conducting channel forms if the electrons near the Fermi 
energy acquire a Berry phase of n due to the translation 
by the Burgers vector when encircling the dislocation. 
This phase twist can be viewed as an Aharonov-Bohm ef- 
fect caused by a fictitious gauge flux $o/2 associated with 
the topological defect and its precise value is protected 
by time- reversal symmetry [9-11]. In two spatial dimen- 
sions, the analogous effect exists: an isolated tt flux in 
a time-reversal invariant TI binds a single Kramer's pair 
which necessarily leads to spin-charge separation [12-15]. 
However, in contrast to three-dimensional systems, the 
realization of this effect by a fictitious gauge-flux gen- 
erated by a topological defect is less universal because 
time-reversal symmetry alone does not protect the pre- 
cise value of the acquired Berry phase [10, 11]. Never- 
theless, spatial symmetries associated with the presence 
of a lattice [16-18] can render the electronic structure of 
a topological point defect robust against local perturba- 
tions also in two dimensions. For example, it has recently 
been demonstrated that dislocations on the square lat- 
tice act sufficiently robust as a &o/2 flux in the topolog- 
ical non-trivial M phase of the Bernevig-Hughes-Zhang 
model for HgTe quantum wells [19] but not in the topo- 
logical T- phase [20] , allowing to distinguish the two topo- 
logically non-trivial phases of this model. 

Motivated by the aforementioned studies, we propose 
in this Letter an alternative way to realize the physics 
of a localized 7r flux in a topological band insulator us- 



ing wedge disclinations on the honeycomb lattice, see 
Fig. 1. Such conical defects form the elementary build- 
ing blocks of various extended lattice defects observed 
in graphene and related carbon based structures [21-26]. 
While the intrinsic spin-orbit coupling in graphene is too 
small to access the TI phase [27, 28] experimentally, sev- 
eral promising routes exist to enhance the spin-orbit cou- 
pling by doping with adatoms [29-31]. We expect that 
the lattice defects discussed in this work are naturally 
present in these graphene-based topological insulators. 
Other potential solid-state realizations of TI phases on 
the honeycomb lattice involve artificial graphene [32] or 
silicene [33, 34]. 

On the hexagonal lattice, an isolated wedge disclina- 
tion is constructed by locally replacing a hexagon by a 
/-polygon (we discuss / = 4, 5, 7, 8) while preserving the 
three-fold connectivity of the honeycomb lattice. This 
introduces a global change of the lattice best illustrated 
by Volterra's cut-and-glue construction [35], in which a 
wedge is removed from or added before gluing the two 
sides back together to form a cone. The point group 
symmetry restricts the possible opening angles to multi- 
ples of 7r/3 and we label different defects by the integer 
n counting the number of removed (n > 0) or added 
(n < 0) 7r/3 wedges. To study the interplay between 
such conical singularities and electrons in topologically 
non-trivial bands, we investigated a model of a Chern 
insulator for spinless fermions on the honeycomb lattice, 
first introduced by Haldane [36]: 

u = -t y £ j (4 c 3 + h - c ) + *a< (''v- r - + h - c ) ■ 

(1) 

The real nearest-neighbor hopping is t and we assume a 
purely imaginary second-neighbor hopping itnV{j where 
Vij = ±1 depends on the hopping direction [36]. At 
half-filling, the model defined in Eq. (1) has a finite 
Hall conductivity a xy = r ye 2 /h with a Chern number 
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FIG. 1. Correspondence between the bound state of (i) and 
isolated -n flux in the defect-free case (n = 0), (ii) an isolated 
pentagon defect (n = 1) with an external flux 7r/2 and (iii) 
an isolated square defect (n = 2) without external flux. The 
external flux $ was applied through the center marked by 
a circle. Results are presented for the model Eq. (f) with 
tult — 0.4. The local density-of-states (LDOS) were obtained 
on a site of the central polygon using the Lanczos algorithm. 



7 = sign(£ 2 i)- A generalization to a time-reversal invari- 
ant topological insulator, in which the imaginary second 
neighbor hopping is generated by intrinsic spin-orbit cou- 
pling, has been discussed by Kane and Mele [27, 28] and 
our results generalize to this situation, as well. 

The main findings of this work is the connection be- 
tween the bound state induced by (i) an isolated $o/2- 
flux, (ii) an isolated pentagon (n = 1) or heptagon 
(n = —1) defect with an external flux of magnitude 
n7$o/4 through the center and (iii) an isolated square 
(n = 2) or octagon (n = —2) defect without external 
flux. We reached this conclusion in three different ways: 
(I) direct computation of the local density of states in 
the lattice model, (II) the analysis of disclinations in the 
continuum model and (III) their description in terms of 
coupled edge modes. 

To compute the local density of states (LDOS) near 



the defect core in the lattice model, we fixed the ra- 
tio tzi/t — 0.4 and used the Lanczos algorithm with 
open boundary condition [37] to obtain the retarded 
local Green's function Gu(E) from which the LDOS 
Ni(E) = -±bD.Gu(E) was derived. We keep up to 300 
states, and use a small imaginary part of 0.02 t to ob- 
tain a smooth spectra. Figure l(i) shows the LDOS in 
the defect-free case on the hexagon through which an 
external flux is threaded. Turning on a finite flux [38] 
moves a bound state from the valence to the conduction 
band reaching E = for $ = $o/2 = 7T. From particle- 
hole symmetry, Ni(E) — Ni(—E), and the conservation 
of states, J_ oo dENi(E) = 1, it follows that the excess 
or deficit charge bound to the ir flux is ±e/2 [12]. Our 
numerical integration of the LDOS confirmed this expec- 
tation. 

Studying the LDOS on a pentagon defect, Fig. 1(h), 
we identify an in-gap state even without external flux. 
Threading $ = it/ 2 through the pentagon shifts the 
bound state energy producing a mid-gap state in anal- 
ogy to the situation (i) with ir flux. On the other hand, a 
flux $ = — 7r/2 produces two symmetric resonances close 
to the band edges. Switching from the pentagon to the 
heptagon defect (n = —1) or changing the sign of the 
Chern number, we find that the opposite sign of the flux 
is required to produce the mid-gap state. 

Figure l(iii) illustrates the case of a square defect. The 
LDOS shows that the mid-gap state is now realized in the 
absence of any external flux. To completely remove the 
bound state, an external 7r-flux is required. We find the 
same behavior also for the octagon defect (not shown). 
Moreover, this property does not rely on particle-hole 
symmetry (or the fact that the bound-state energy is at 
E = 0): in an analogous calculation including also real 
second-neighbor hopping, we find that only an external 
Tr-flux is able to completely remove the bound state. The 
robustness of the correspondence between (i), (ii) and 
(iii) is further discussed below. 

The numerical results presented in Fig. 1 can be con- 
sistently explained from a continuum description, as we 
discuss in the following. In the low-energy limit, Eq. (1) 
reduces to the Dirac Hamiltonian with a "Haldane mass" 
m = 3\/3i 2 j: 



H = v[r z a x p x + dyPy] + mr z a z 



(2) 



where a = (u x ,a y ,a z ), f = {t x ,t v ,t z ) are Pauli ma- 
trices denoting the sublattice and valley degrees of free- 
dom, respectively, and v = y/3ta/(2h) is the Fermi ve- 
locity. H acts on the four-component spinor ^(r) = 
[ipA(r), %[)B(r), ifjA> (r), ipB'(r)] T where A and B label the 
sublattice in valley K and A' and B' in valley K' [38] . It 
is known that deformations of the honeycomb lattice en- 
ter the continuum description via fictitious gauge fields 
[39-41]. As we review below, topological point defects 
manifest themselves by spatially well-localized fluxes of 
the fictitious fields [21]. 
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FIG. 2. (a) Continuum version of the cut-and-glue construc- 
tion with a regularization hole of radius p around the origin. 
a(8) is a closed path around the cone, (b) The boundary 
conditions for the spinor across the seam have to compensate 
the mismatch of the base functions, as indicated for a pair 
of matching degrees of freedom A and B' for the pentagon 
disclination. 



We model the disclination by a regularized cone where 
a disk of radius p around the apex is removed, see 
Fig. 2(a). The fictitious gauge fields are related to the 
non-trivial holonomy when the spinor is parallel trans- 
ported along a closed path a{9) (0 < 9 < 2ir) around the 
cone: 



*(a(27r)) =ft n tf(a(0)), U n = e 



0. 



(3) 



Tin simply represents the rotation by the Frank angle 
mr/3 of the defect. The boundary condition for the en- 
velope function, Eq. (3), compensates the mismatch of 
the base functions e ±lK r across the seam in the cut- 
and-glue procedure illustrated in Fig. 2(b), making the 
total wave function single- valued [21, 22, 38, 39, 42]. 

To deal with Eq. (3) for general n, we seek for a local 
gauge in which the boundary condition is independent of 
disclination type. To achieve this goal, we introduce po- 
lar coordinates (r, <fi) defined in the unfolded plane, see 
Fig. 2(a), and perform two singular gauge transforma- 



tions * H> * i-4 with 



V n {6) = e iS z- a « T v. 



(4) 



The first operation U transforms <P to a co-rotating 
spinor [43], effectively replacing d r by d r + l/(2r) in the 
Hamiltonian. The second gauge transform V n introduces 
a matrix-valued gauge field into the Hamiltonian, effec- 
tively replacing dg by dg — ijO- y T y . The transformed 
spinor ^> n (r,9), 9 = 4>/(l — |), is now anti-periodic in 
9 for any n. In the final step, we use a global transfor- 
mation S 



s = 



1 

71 



(l + iT x (Ty), (5) 



which block-diagonalizes the Hamiltonian and defines 
two emergent valleys t = ±. The separation into two 
decoupled valleys is not always possible in the presence 
of a mass term but reflects the fact that the Haldane 
mass m,T z a z is compatible with all types of disclinations 
[44] . The block-diagonal Hamiltonian H' n = UHW with 
U = SV n U, has the same form for any n. The product 



ansatz &' n (r, ff) = x(r)e 4:,e with half-integer j decouples 
radial and angular part and the radial part is (Jl = 1 = v) 

[381 



H' T (n) 



rd r + i ) to x + w T {n)o y 



mra z 



(6) 



As before, r = ± denotes the emergent valley and [22, 23] 



y T (n) 



1 - § 



(7) 



Equation (7) also accounts for a localized real magnetic 
flux $ through the origin [38]. The topological defect 
manifests itself through the denominator 1 — n/6 and an 
additional gauge flux of 7ir$ / 4 with opposite sign in the 
two emergent valleys. 

The eigenvalue problem H' t (ji)xt{t) = Exrir) can be 
solved in each valley separately. Bound states are given 
in terms of modified Bessel functions of the second kind 
which decay exponentially for r — > oo. We find 



X+(r) = 
X-(r) = 



K{ v+ -l/2){Kr) 
-i^K {v _ +1/2) {KT) 



(8) 
(9) 



where K = v m 2 — E 2 > 0. The square integrability of 
^ for p — » does not uniquely determine the bound 
state. To obtain quantized solutions, the internal struc- 
ture of the disk r < p has to be specified. The correct 
quantization is achieved by replacing the Haldane mass 
in Eq. (6) by a confining potential V(r < p) = —Ma z 
[45, 46]. The mass term — Ma z , as compared to the Hal- 
dane mass, has opposite sign in one of the emergent val- 
leys, thereby defining the topologically trivial insulator. 
Because U^aJJ — o- x t x , we identify V(r < p) in the 
frame of Eq. (2) with the inversion symmetric mass term 
of a kekule distortion [47, 48]. In the limit M — > +oo, 
matching of the wave function at r — p takes the form of 
a linear restriction [49] 



7^-j)*(p > ^)=*(p > 0). 



(10) 



r y ) T is the normalized axial current 



Here, J = (— t x o 
in the frame of Eq. (2) and the azimuthal unit vec- 
tor. The sign on the left-hand-side is fixed by the Chern 
number 7 = sign(m). As expected, for the spinors of the 
form Eq. (9), the infinite-mass condition Eq. (10) can be 
satisfied only in the valley with the sign change of the 
mass. Moreover, Eq. (10) leads to the quantization of 
the bound-state energy through 



m 



E K^ y _ 1/2 (Kp) 



m 



E K 



v-, + l/2(Kp) 



(11) 



which, in combination with Eq. (7), incorporates the 
main result of the present work. As illustrated in 



4 



Fig. 3(a) for different values of the dimensionless ra- 
dius p\m\, Eq. (11) has monotonic real solutions for the 
bound-state energy \E\ < \m\ as function of v in a range 
\v\ < 1/2 + p\rn\. For u = it follows that E = inde- 
pendent of p\m\. The physical bound-state spectrum as 
function of flux $/<E>o for a specific defect is constructed 
from the general solution by use of Eq. (7) and is found 
to agree with the numerical results obtained in the lat- 
tice model. The case of a pentagon defect (n = 1) with 
either sign of the Chern number 7 = ±1 is illustrate in 
Fig. 3(b). In particular, this solution predicts that in- 
sertion of an external magnetic flux 771-/2 [marked with 
o] shifts the bound state to zero energy while the op- 
posite flux — 77r/2 [marked with □] leads to two bound 
states symmetrically arranged with respect to E = 0, in 
accordance with the results shown in Fig. l(ii). 
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FIG. 3. (a) General solution for the bound-state energy from 
Eq. (11) as function of v for different radii of the hole p\m\. 
(b) Bound-state energies as function of external flux $/<&o for 
the pentagon defect. 

The correspondence between external magnetic and in- 
ternal fictitious fluxes induced by wedge disclinations has 
an intuitive explanation via the coupling of edge modes 
across the seam. This picture also gives a deeper under- 
standing of the robustness of our results. Our analysis 
follows related work for dislocations [9, 11]. We discuss 
the square defect but the argument applies equally-well 
for the other wedge disclinations. Let us imagine two 
disconnected flat honeycomb sheets from which a 120°- 
wedge has been removed, as illustrated in Fig. 4(a). The 
bulk-edge correspondence for Chern insulators implies 
chiral edge states propagating along the zig-zag edges 
of top and bottom part. In the vicinity of the energy 
crossing, they are described by the edge theory 



Hedge = J d£<pt(£)[-ivdt<T, + (12) 

with /i(£) = 0. The wave function <ys(£) varies smoothly 
on the scale of the lattice constant and includes right and 
left moving components p = (pn,tpL) T ; £ is the coordi- 
nate along the cut. The total edge wave function on a lat- 



tice site is given by v 5 cd g e(£) 



-ik B £ 




(b) EA 



2tv 
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where kg is the edge momentum. Inversion symmetry 
implies that the edge states of a zig-zag edge cross at 



FIG. 4. (a) Coupling of chiral edge states across the seam in 
the cut and glue construction of a n = 2 disclination. (b) In 
the absence of a coupling between top and bottom part, the 
edge states cross at fed = ir, which is protected by inversion 
symmetry, (c) Turning on a weak coupling introduces a in- 
dependent mass term /x(£) where £ is the coordinate along 
the edge. This opens a gap in the edge spectrum (d). The 
modified matching condition for £ < and £ > effectively 
leads to a sign change of the mass /x(£) which hosts a bound 
state. 



either fcga = or ir - in the model Eq. (1), they cross 
at kEd = 7T, see Fig. 4(b). Hence, the base functions 
giafcjsC oscillates with a period of two, as indicated in 
Fig. 4(a). Imagine now to weakly couple left and right 
movers which is described by p ^ in Eq. (12). This 
gluing across the seam opens a gap of order p in the edge 
spectrum, Fig. 4(d). However, to account for the differ- 
ence between £ > and £ < in the matching of the base 
functions, p(£) has to change sign at £ = 0, Fig. 4(c). It 
is this sign change which is responsible for the existence 
of a bound state, in analogy to solitons in polyacetylene 
[50]. The sign change of p is equivalent to introducing a 
ir flux in the defect- free system [12]. 

If inversion symmetry is broken, the edge states cross 
away from fc^a = ir. However, as long as the edge state 
theory can be obtained by expansion with the base func- 
tions at kEa — 7T, the correspondence to an isolated ir flux 
remains, even though the bound-state energy in general 
shifts away from E = 0. We have numerically confirmed 
this expected robustness by adding both local perturba- 
tions in the form of on-site potentials as well as various 
global symmetry-breaking terms including staggered sub- 
lattice potentials, real second neighbor hopping as well 
as dimerized first-neighbor hopping. 

Our main results are summarized in Fig. 1 and given 
by Eq. (11) in combination with Eq. (7) which establish 
a correspondence between an external magnetic flux and 
internal fictitious fluxes of topological defects in a Chern 
insulator on the honeycomb lattice. Our results also gen- 
eralize to time-reversal invariant topological insulators. 
The disclinations then act as a source of spin-flux [14], 
i.e. a flux with opposite sign for the two spin components. 
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We suggest that the characteristic spectrum induced by 
the topological defects can be measured by scanning tun- 
neling microscopy. 
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FIG. SI. The complex phases of the base functions in graphene. A and B refer to the sublattice components at valley K while 
A' and B' refer to the sublattice components at K' . We use the notation rj = e 27rl//3 and fj = e -2 ^^ 3 . (a) shows u A {r) and 
u B i(r) while (b) shows u A /(r) and ub(t). 



Supplementary materials 



I. BASE FUNCTIONS AND BOUNDARY CONDITIONS 



The boundary conditions in the cut and glue procedure for the different disclinations are obtained by matching the 
total wave function across the cut [22]. To derive the continuum theory, this requires to fix a convention for the base 
function. It is convenient to choose a set of base functions which explicitly preserves the inversion symmetry with 
respect to the origin located at a center of a hexagon. For valley K = (An/ (3V3d), 0), the base functions are denoted 
by u A (r) and ug(r); for valley K' = — K they are denoted by u A /(r) and u B >(r). The complex phases are given by 
the Bloch factors e ±lr ' K with r the location of the given A or B site. Using the notation rj = g 2, ™/3 anc [ fj — e -2iri/3^ 
the amplitude of the base functions are given in Fig. SI. The total wave function is then expanded around the four 
base functions as 

Aot{r) = ip A (r)u A (r) +ip B (r)u B (r) + i(> A >(r)u A >(r) + ip B >(r)u B ,(r), (SI) 

where we have suppressed the dependence on the physical spin. The envelope function 

*(r) = [rf>A{r),i> B (r),Mr),i>B'(r)] T (S2) 

is a four-component function which is assumed to be smooth on the atomic scale. In this convention, the Dirac 
Hamiltonian in the presence of a Haldane mass is given by Eq. (2) 

H = v(r z a x p x + a v py) + mr z (j z (S3) 

where a = {cr x ,cry,a z ) and f = (T x ,T yi T z ) are Pauli matrices denoting the sublattice and valley degrees of freedom, 
respectively. 

The boundary conditions of the envelope function are to make sure that the phase miss-matches of the base functions, 
as shown in Fig. S2 and S3, are properly compensated so that the total wave function is single- valued upon encircling 
the defect core. For the even-membered disclinations (n — ±2), the single valucdncss of the total wave function 
implies the following boundary conditions for the envelope function: 

^(r,<p = 4:Tr/3) = e- i ^ T ^(r,(t> = 0), n = 2, (S4) 

*(r,<£ = 8tt/3) = e +l ^ a * T ^{r,<j) = 0), n = -2. (S5) 

Note that these boundary conditions are diagonal in sublattice and valley degrees of freedom. 

For the ±60°-disclinations the sublattice components are no longer conserved. Nevertheless, it is possible to find 
boundary conditions which are independent of the distance from the defect core by matching opposite sublattice 
components in opposite valleys, as shown in Fig. S3. The boundary conditions are then given by 



(0 

fj 

fj 

\t} 0/ 



*(0 = O), n = l, (S6) 



FIG. S2. The amplitudes of the base functions of valley K across the cut for the 120°-disclination. The amplitudes for valley 
K' are obtained by complex conjugation of the above shown amplitudes. 




FIG. S3. The amplitudes of the base functions corresponding to the components which are matched across the cut for the 
60°-disclination. The remaining two components are obtained by complex conjugation of the above shown amplitudes. 



for the pentagon (+60°) disclination and by 



7tt 



fO fj\ 

?y 

rj 

\fj 0/ 



$(0 = 0), n = -l, 



(S7) 



for the heptagon (—60°) disclination. With the rescaled angular variable 8 — <fi/£l n (f2„ = 1 — n/6), Eq. (S5)-(S8) can 
be written in a compact form as Eq. (3): 



2tt) = e 3 



'*(0 = 0). 



(S8) 



II. BLOCK DIAGONALIZATION OF DIRAC HAMILTONIAN 



After performing the transformation to the co-rotating frame [see Eq. (4)], 

( H Q i K ) = U W HUt W> UW = e i t<* T >, (S9) 
the Dirac Hamiltonian in polar coordinates reads 



id r — -d, 



■<t> 



f2„*n 2 



( + - 

\n n <S>„ ~ 2 



H K =\ , N r r 1 . (S10) 



for valley K. In the opposite valley, it is given by 



n I ~ m ^-F^-Hra~i| (SII) 



We have also included a magnetic flux $ through the origin, see next section. We introduce the rescaled angular 
variable = <p/Q n (f2 n = 1 — n/6), and perform the second local transformation with 



V„(H) = < ayTy = cos ( ^j-J cr T + i(a v T v )sm 



H n = V n (8)\ H « ~° )v}(0). 



to obtain 

H„ = V„(0\ ( J 

H K > 

After these two gauge transformations, the transformed spinor 

*„ = V n U9 

now obeys anti-periodic boundary conditions for any n, 

* n (0 = 2tt) = -9 n {6 = 0), 
instead of the awkward condition Eq. (3). The transformed Hamiltonian explicitly reads 



(S12) 



(S13) 



(S14) 



(S15) 



H n = 



— z n 
rO n 4 



— m 


— i n 
rfl„ 4 



2 n 
rf2„ 4 



rO„ 4 



It is brought to block-diagonal form by applying the global transformation 

SHnS^, S = -j= (1 + ir^o-y) 

Separation of angular and radial variables, 

<M)= X (r)e^ 
with j half-integer to satisfy Eq. (S15), leads to Eq. (6). 



(S16) 
(S17) 

(S18) 



III. MAGNETIC FLUX WITH A DISCLINATION 



Continuum model 



In the effective Hamiltonian, a magnetic flux <I> passing through origin is described by replacing the momentum p 
in Eq. (2) by the canonical momentum p + A with the vector potential 



A(r,< 



$ 1 



■(— sin^, cos< 



(S19) 



2-7T f2„r 

with Q„ = 1 — n/6. In polar coordinates, the magnetic flux enters the Hamiltonian as shown in Eqs. (S10) and (Sll 



Lattice model 



In the lattice model, a magnetic flux $ through the origin is described by modifying the complex phase of the 
hoppings according to the Peierls substitution. The choice of complex phases is not unique (since many vector 
potentials lead to the same magnetic field), and one simple arrangement is now provided: one draws an arbitrary 



9 



semi-infinite string starting from the origin, and all hopping intersecting with this string is attached with the phase 
<3> or -$ (e = h = 1) via 

Uj -> ^. e ** si s n [( R *- R ,) x M (S20) 

where ijj (generally a complex number) is the hopping amplitude from site i at to site j at Rj in the absence 
of external fluxes, and t is the tangent of the string at the intersection. In our numerical simulations we choose the 
semi-infinite string to be a straight line starting from the disclination center and crossing the middle of one of the 
edges of the core polygon. 



